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We present a novel technique for the calculation of dynamical correlation functions of quantum 
impurity systems in equilibrium with Wilson's numerical renormalization group. Our formulation 
is based on a complete basis set of the Wilson chain. In contrast to all previous methods, it 
does not suffer from overcounting of excitation. By construction, it always fulfills sum rules for 
spectral functions. Furthermore, it accurately reproduces local thermodynamic expectation values, 
such as occupancy and magnetization, obtained directly from the numerical renormalization group 
calculations. 
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I. INTRODUCTION 

The understanding of quantum impurity systems be- 
came one the cornerstones of condensed matter theory 
during the last decade. Quantum impurity systems ap- 
pear at the heart of a variety of different physical prob- 
lems. Traditionally, they were used to describe the in- 
teraction of magnetic impurities with a metallic host. 1 
Nowadays, quantum impurity systems play a fundamen- 
tal role in our understanding of the low temperature 
properties of single-electron transistors^ and the tun- 
neling spectroscopy of adatoms on metal surfaces4£ The 
basic structure common to all such systems is a meso- 
scopic subsystem or device - such as a quantum dot, an 
organic donor-acceptor molecule or an adatom - coupled 
to a continuum of states that can be represented by non- 
interacting particles of cither fcrmionic or bosonic na- 
ture. Typical realizations are models like the single im- 
purity Anderson model^ where the continuum of states 
is described by free fermions, or the spin-boson model, 7 
where the discrete orbitals interact with a bosonic bath. 
In addition, within the dynamical mean-held theory&S 
or its cluster extensions^ lattice models for strongly cor- 
related fermions have been mapped onto quantum im- 
purity problems embedded in a fictitious, self-consistent 
bath. Obtaining a self-consistent solution for these theo- 
ries requires very accurate determination of local Green's 
functions. 

Wilson's numerical renormalization group (NRG) 
approachii to quantum impurity problems is one of the 
most powerful and flexible ways for accurately calculat- 
ing thermodynamic properties of quantum impurity sys- 
tems. In addition, it provides a deep insight into the un- 
derlying physics through the analysis of the fixed point 
Hamiltoniansii2iiiii4 The key ingredient is a logarithmic 
discretization of the bath continuum, resulting in a well 
defined hierarchy of energy or temperature scales. The 
discretized model is then iteratively diagonalized and the 
basis set truncated, retaining only those states with low 



lying energies after each step. Each iteration represents 
a certain temperature T, and all thermodynamic proper- 
ties are determined for that particular Tiii*i£ For calcu- 
lating dynamical propertiesi^iSiiLi2iiS three fundamen- 
tal problems arise: (i) how to recover the continuum limit 
from a discretized spectrum, (ii) how to obtain dynamical 
information on all energies scale at arbitrarily low tem- 
perature in such away that (iii) spectral sum rules are 
always fulfilled and thermodynamics expectation values 
are reproduced exactly independent of how many states 
are kept after each NRG step. 

In this paper, we will derive an exact analytical expres- 
sion for arbitrary dynamical correlation functions which 
solves this problems. It is based upon the recent iden- 
tification of a complete basis set of the Wilson chain, 
which is also an approximate eigenbasis of the NRG 
Hamiltoniani22i2i We will show that this complete ba- 
sis set automatically ensures that spectral sum rules are 
fulfilled exactly and thermodynamic expectation values 
reproduced accurately by the spectral functions. We 
furthermore explicitly demonstrate that our novel ap- 
proach yields spectral functions which are largely insen- 
sitive to the number of eigenstates kept in each NRG 
iteration, thus improving the applicability of the NRG to 
multi-orbital and multi-site impurity problems and con- 
sequently also to multi-band and cluster DMFT prob- 
lems, because the number of NRG states needed can be 
significantly reduced without loosing accuracy. In ad- 
dition, phenomenological patching algorithms typically 
employed to combine excitations from different NRG 
iterations^ become obsolete, since we rigorously identify 
which excitations actually contribute at which Wilson 
shell. 

The paper is organized as follows. In the next section 
we will briefly review the basic theoretical concepts and 
derive the complete basis set. We furthermore show how 
this complete basis set can be used to calculate dynamical 
correlation functions. Section III is devoted to a detailed 
discussion of results for the simplest and most important 
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quantum impurity model, the single impurity Anderson 
model. We will explicitly compare our new method to 
the standard implementations. A summary and outlook 
in section IV will conclude the paper. Proofs for our 
claim that sum rules and spectral averages are automat- 
ically respected within this formulation are given in the 
appendix. 

II. THEORY 

The NRG is a very powerful tool for accurately cal- 
culating equilibrium properties of quantum impurity 
models. Originally developed for treating the single- 
channel, single-impurity Kondo Hamiltonian^ 1 * 2 * this 
non-perturbative approach was successfully extended to 
the Anderson impurity model^^S^ the two-channel 
Anderson^ 3 - and Kondo Hamiltonian, 24-25 different two- 
impurity clusters j 26 i 27 i 28 i 29 i 30 i 31 i 32 and a host of related 
zero-dimensional problems. Recently, it was extended to 
equilibrium properties of impurity models with a bosonic 
bath 3 ^ 3 ! or even combinations of both fermionic and 
Bosonic baths & 

The Hamiltonian of a quantum impurity system is gen- 
erally given by 

H = 7ibath + Wimp + Wmix j (1) 

where 7ibath models the continuous bath, Hi mp represents 
the decoupled impurity, and H mlx describes the coupling 
between the two subsystems. Thermodynamic properties 
of such a quantum impurity system are very accurately 
obtained using the NRG. At the heart of this approach is 
a logarithmic discretization of the continuous bath, con- 
trolled by the discretization parameter A > ljil the con- 
tinuum limit is recovered for A — > 1. Using an appropri- 
ate unitary transformation^ the Hamiltonian is mapped 
onto a semi-infinite chain, with the impurity coupled to 
the open end. By construction, the A th site of this chain 
couples only to its immediate neighbors, which allows to 
write the Hamiltonian of the infinite chain as limit for 
N — > oo of a sequence of finite chains, denoted by Hm, 
with a unique prescription Hm i— * Hm+i, the RG equa- 
tion of the NRG. The ./V th link along the chain represents 
an exponentially decreasing energy scale D m ~ Ar N / 2 for 
a fermionic 4 i and Dm ~ A~ N for a bosonic bathi^ 3 . Using 
this hierarchy of scales, the sequence of finite-size Hamil- 
tonians Hm for the A-site chain 4 ! is solved iteratively, 
i.e. starting with N = one constructs the Hamiltonian 
Hq, diagonalizes it, adds the next site to obtain H±, diag- 
onalizes it etc. The problem of an exponentially increas- 
ing Hilbert space is resolved by observing, that, due to 
the exponential decrease of the energy scales, only low- 
energy states actually contribute in the step N — ► N + 1 ; 
one thus discards the high-energy states before moving 
on to the next step to maintain a manageable number of 
states, we denote with Ns in the following. The reduced 
basis set of Hm obtained that way is expected to faith- 
fully describe the spectrum of the full Hamiltonian on a 



scale of the order of Dm, corresponding to a temperature 
Tjv ~ Djviii Note that even from this brief discussion it 
is evident that for a given step N all information about 
energies E ^S> Dm has been lost, while no information 
about energy scales E <C Dm is available yet. Thus, to 
calculate dynamical quantities with a similar accuracy 
as thermodynamic, one has to tackle the problem of cor- 
rectly mixing information from earlier and later NRG 
steps. 



A. Complete Basis Set 

In a recent extension of the NRG to real-time dynam- 
ics out of equilibriumSSiSi a complete basis set for such 
a Wilson chain of length N has been identified. Further- 
more, this complete basis set also forms an approximate 
eigenbasis of the Hamiltonian Hm- Since this complete 
basis set plays a crucial role in the derivation of the an- 
alytical expression for spectral functions, we summarize 
briefly the main ideas^S of the proof of completeness dis- 
cussed extensively in Ref. |2l[ 

At first sight the claim that the NRG automatically 
generates a complete basis set, which is also an approx- 
imate eigenbasis of the chain Hamiltonian Hm, might 
appear contradictory. A renormalization procedure is 
usually viewed as a clever way to identify the relevant 
degrees of freedom by reducing the dimensionally of the 
underlying Fock space. In order to excavate the optimally 
adapted complete basis set, we have to shift perspectives 
how to view the NRG algorithm. 

There are two possible ways to interpret the iterative 
NRG solution of the TV-site chain. In the traditional pic- 
ture one starts from a core cluster that consists of the 
impurity degrees of freedom and the A = site, and 
enlarges the chain by one site at each NRG step. Al- 
ternatively, one can view the NRG procedure as starting 
from the full chain of length N, but with the hopping 
matrix elements set to zero along the chain. At each suc- 
cessive step another hopping matrix element is switched 
on, until the full Hamiltonian TLm is recovered. In this 
latter picture, to be adopted below, the entire sequence 
of Hamiltonians H m with m < N act on the Fock space 
of the A-site chain. 

Accordingly, each NRG eigen-energy of H m has an ex- 
tra degeneracy of S N ~ m \ where d is the number of dis- 
tinct states at each site along the chain. The extra de- 
generacy stems from the N — m "environment" sites at 
the end of the chain, depicted in Fig. ^ 

The set of eigenstates of H m , conventionally denoted 
as {|r)}, can be formally constructed from the complete 
basis set {|ai mp , cxq, ■ ■ ■ , a at)} of the NRG chain of length 
N where the on label the configurations on each chain link 
i. Since H m does not act on the chain links m+1, • • • , A, 
\r) can be written as \r,e;m) where the "environment" 
variable e = {a m+ i, ■ ■ ■ ,o:n} encodes the A — to site la- 
bels a m +i, • • • ,olm- The index m is used in this notation 
to record where the chain is partitioned into a "subsys- 
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FIG. 1: (color online) The full Wilson chain of length N is 
divided into a sub-chain of length m and the "environment" 
Rm,N- The Hamiltonian TL m can be viewed either as act- 
ing only on the sub-chain of length m, or as acting on the 
full chain of length N, but with the hopping matrix elements 
t m , • • • ,ijv-i all set to zero. The former picture is the tradi- 
tional one. In this paper we adopt the latter point of view. 



with the completeness relation 

1 = l m + l£ • (5) 
Note that for m — N only 1~ exists. 

B. Definition of the Green's function 



tern" and an "environment" (see Fig.^). 

Consider now the first iteration m m ; n at which states 
are discarded. In order to keep track of the complete 
basis set of the iV-site chain, we formally divide the 
eigenstates \r,e;m) into two distinct subsets: the dis- 
carded high-energy states {\l, e; m)dis} and the retained 
low-energy states {\k,e;m)k P }- In the course of the 
NRG, only the latter states are used to construct the 
next Hamilton matrix Hm+i within the reduced subspace 
{|fc, a m +i, e'; m)}. Note, however, that the sum of both 
subsets still constitutes a complete basis set for the N- 
site chain. Repeating this procedure at each subsequent 
iteration, we recursively divide the retained subset into 
a discarded part and a retained part. Then, the collec- 
tion of all discarded eigenstates \l, e; m)di S together with 
the eigenstates of the final NRG iteration N combine to 
form a complete basis set for the entire Fock space Tn- 
Regarding all eigenstates of the final NRG iteration as 
discarded, one can formally write the Fock space of the 
iV-site chain in the form Tn — span{ 1 1, e; m)dis}- Since 
all states are retained in the course of this construct, the 
following completeness relation obviously holds: 



JV 

E 



\l,e;m)dis dis(l,e;m\ 



(2) 



Here the summation over m starts from the first iteration 
"imin at which a basis-set reduction is imposed. All traces 
below will be carried out with respect to this basis set. 
Hence, the evaluation of the spectral functions will not 
involve any truncation error. Note also that we made no 
reference to a particular Hamiltonian H in constructing 
the basis set {|Z, e; m)dis\- 

Another useful identity to be used below pertains to 
the subspace retained at iteration m, {\k, e; m}k P }- To 
this end, we note that the sum in Eq. J5J can always be 
divided into two complementary parts 1~ and 1+ : 



i; 



^2\l',e';m') dis dis (l',e';m'\, (3) 



m' —m m i n I' 
N 



E y^J l ',e';m') dis dis (l',e';m'\ . 

m'— m+1 ,e' 

\k, e; m) kp kp (k, e; m\ . (4) 



The textbook definition of the retarded Green's func- 
tion is given by 



G AtB {t) = -iO(t)Tr[p[A(t)B}_ 



(6) 



where [A, B] s = AB — sBA with s = 1 for bosonic op- 
erators A, B and s = — 1 for fermionic operators. In our 
investigation, we restrict ourselves to local operators A 
and BM 

As already mentioned, the fundamental philosophy of 
the NRG is that a chain of length N corresponds to a 
temperature scale T ~ D^, hence [3E" 1 3> 1 for all m < 
N, where Ef 1 denotes the eigen energies of H m - In the 
NRG, the thermodynamic density operator p is therefore 
represented only by the states of the last iteration N 



P = 



Pi 



Z 

Zn 



^ H ^^-J2 Pl \l;N){l;N\ (7) 

(8) 



and Zn — J2i e 

If we were able to solve the NRG chain of length N 
without any truncation, the Green's function Ga,b(z) 
would be given by the textbook Lchmann representation 



Ga,b(z) 



/ dte lzt G A . B {t) 
Jo 



z + Ei- E v 



(9) 



Derivation of the NRG Green's function 



Let us start from the retarded Green's function looking 
only at the first term from the commutator and insert two 
completeness relations (0. We then obtain 
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Tr [pe im Ae- im B] = ^ ^ Tr [pe lHt \l, e; m)(l, e; m\Ae~ im \l\ e'; m')(l', e'; m'\B] 

l,e,m V ,e' ,m' 

= e; ™l^ e " ifft l z '> e '; diaC'i e '; m|-B/5e ifft |/, e; m) dis 

m Z,e C,e' 

+ X! X! * s ^' e; TO l^ e_lfft e '; TO )fep fcp( fc > e '; m\Bpe lHt \l, e; m) rfis 

m Z,e k.e' 

+ dls (l,e';m\Be lHt p\k,e;m) kpkp (k,e;rn\Ae~ lHt \l,e';m) dls . 

m J,e' k,e 



For the second term from the commutator in JJJJ one 
obtains a similar expression. Here and in the following 
we use the convention that an index I labels a discarded 
state, while an index k represents a state kept at a certain 
chain length m. The first term in (|10f> arises from m! — 
to, the second from to' > m and the third from ml < 
to, where we made use of Eq. (0J. This trick 20 allows 
to express the Green's function as sum over equal shell 
contributions only. Note that this is an exact formula: 
no approximations have been made so far! 

Since the state \s,e;m) is an eigenstate of H m , i.e. 
H m \s,e;m) — E™\s, e;m), we will now use the approx- 
imation H\s,e;m) w E™\s, e;m). This approximation, 
justified due to the energy hierarchy implied by the log- 
arithmic discretization^ is completely in the spirit of 
the NRG and in fact used in the calculation of ther- 
modynamics properties*^ Note, that this will be the 
only approximation made, which is of energetic nature 
and unrelated to any truncation error. We thus eval- 



(10) 

I 

uate exp(iHt)\s, e;m) ~ exp(iE™t)\s, e; m) and Laplace 
transform all contributions. The first term from (|10|) and 
the corresponding expression for the second term from 
the commutator contains only discarded states and re- 
duces, therefore, to the last iteration m = N due to the 
representation of the density operator in NRG by Eq. J7J). 
It yields as contribution to Ga,b{z) 

X z + E^-EP 

In the next two terms the summation over all en- 
ergy shells has to be evaluated, because the states 
\k, e; m)k p are not orthogonal to \l\ N). We obtain, using 
p\l, e; m) dls = for to < 



Ga,b( z ) 



N-l 



m—mmin l,e k 



> 2 Mi e ' m \Ap\k, e';m)(k,e';m\B\l,e; 



—s 



z + E™- E™ 



(12) 



and 



N-l 



G% B (z) = J2 EE( , ' e 'i m l%ei'»HM;HA|i 1 e'; 



m=m mi „ Z,e' k,e 



z + E™~ E" 



(13) 



Since in the last iteration there are no kept states, all 
are considered discarded, to = N does not contribute 
to G" and G m . Now we insert a completeness relation 
(1+ +1~) between A and the density operator in G % \ B (z) 
and B and the density operator in G™ B (z) and make use 
of l~p = for all to < iV. 21 Because A and B are local 
operators^ the matrix elements 



are diagonal and independent of the environment vari- 
ables e. n e > denotes the number of Fermions in the en- 
vironment times the total number of Fermions created 
by A. Since the matrix elements of A and B are eval- 
uated simultaneously, the total phase factor is given by 
[s n ''} 2 — 1 for operator A and B^ imposing the same 
change of the total particle number. Therefore, the trace 
over the environment only acts on the density operator 



(k, e; m\A\l, e'; to) = s n -' d e ^A kd (m) (14) 
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p and the final result lAlthat the spectral sum rule 

N - 1 r dz 

G% B {z) = Yl EE A '."'( m )^( m ) 5 wH C = f^n GAM 



l k,k' 

-S 



z + E[-E k 



(15) 



£ f^-G%B{*) = 'to\p[A,B].]) (19) 



a— 



an< ^ is always fulfilled exactly. 

jv-i Therefore, the spectral function will be correctly nor- 

G% B {z) = V V y2B lyk ,{m)pl?, d k {m)A k j{m) malized independent of the number of states N s kept 

~ ' L ' ' after each NRG iteration. Note that the NRG trunca- 



n l k.k' 

1 



tion only influences the partitioning of the states, but 
(16) never the completeness of the basis. Hence, the spectral 
z k ... functions become more robust to truncation errors, as 

Ga,b(z) — G l A B (z) + G^ B (z) + G^\b( z ) (17) will be shown in section ITTT1 by using an extremely low 

number of kept NRG states. The spectral sum rule is 
is formulated in terms of the reduced density matrbi2£i21 a consequence of an operator identity and independent 

of approximations made in the dynamics as long as no 
states arc discarded in the calculation. 



PM'M^E^H^e;™) , (18) 



first introduced to the calculations of NRG spectral func- 
tions by Hofstetterpi^ The part G"(z) describes negative 
frequency excitations, while G m (z) accounts for positive 
frequency excitations for all m < N because Ei — E k > 
by construction. G l (z) sums all excitations of the last 
iteration iV and has the form of the usual Lehmann rep- 
resentation Eq. ©. 

A few words are in place to pinpoint the difference to 
the Hofstetter approach-^ to calculate spectral functions 
with NRG. Although our reduced density matrix p rcd is 
identical to the one given in Eq. (7) in Ref. Q our rig- 
orous derivation differs in the summation of excitations 
contributing to the Green's function. Eqs. i|15|) and l|16l) 
state clearly that one must only include excitations be- 
tween a discarded and kept state while in the Hofstetter 
approach 1 — the summation index I runs over all states 
present at iteration m. This leads to an overcounting 
of contributions, by the way inherent to all previous ap- 
proaches to NRG spectral functions, see for details Refs. 
ll5Bl6lll7lll8lll9l and references therein. The origin of the 
restriction of summation in Eqs. (|15f) and l|15|) is quite ob- 
vious: All kept states in iteration m span the Fock-space 
of Hamiltonian H m+ i and, therefore, will contribute to 
the excitations at a latter iteration m! > m. They must 
not be included at iteration m. 

In addition, the phenomenological patching 
algorithms, 1 — where spectral information from dif- 
ferent energy shells are "merged" become obsolete in 
our approach. Equations (|1 111 tj|> state exactly which 
excitations contribute at which Wilson shell. Obeying 
this summation restriction ensures the fulfillment of the 
spectral sum rules independent of the number Ns of 
NRG states kept at each iteration. With a little algebra 
and the completeness relations (JSJ) , we show in appendix 



D. Occupation number 

Let us specialize to the local spin-dependent fcrmionic 
spectral function, i.e. A — f a and B — /J. The local 
occupation (/J/o-) can be expressed with the spectral in- 
tegral as 

/OO 7 
— f(u)SmG uj} (u-i5) 

= ftl I{z)G f^ z) (20) 

where f(w) is the Fermi function. By substituting Eqs. 
(I11I16II for G f f \(z) and evaluating the contour inte- 
gration, we show in appendix |B] explicitly that our ex- 
pressions approximately - at T = even exactly - re- 
produce the expectation value (flfo) calculated directly 
with the NRG. In contrast to the spectral sum rule, how- 
ever, which is exact and, therefore, reproduced in our 
numerics with machine accuracy, the accuracy of the oc- 
cupation numbers calculated from the spectra depends 
on the validity of the assumption of vanishing Boltzmann 
factors. Therefore, we expect it to show a certain error at 
finite temperatures. Nevertheless, we find that the devi- 
ation between the NRG values and the ones obtained by 
the numeric evaluation^ of H2(J|I remains less than 10~ 4 . 
This implies also that the NRG value for quantities like 
the local magnetization m = cr(/ ( |/ cr )/2 is accurately 
reproduced by our formulation of the Green's function. 
This is a significant improvement over the Hofstetter ap- 
proach, where deviations on the percent level between 
m NRG anc j m GF nave been reported. 1 - 
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III. RESULTS 
A. Single impurity Anderson model 

The general scheme for the calculation of spectral func- 
tions presented in section [H] did not make any reference 
to a certain model or bath statistics. Therefore, our al- 
gorithm is suitable for any quantum impurity system. In 
order to demonstrate the virtue of the new approach, we 
will present in this section as an important example cal- 
culations for the spin-dependent single-particle spectral 
functions of the single impurity Anderson model (SIAM) 
with and without an external magnetic held at T = 
and hnite T . 

The Hamiltonian of the SIAM 6 i 12 i 13 

ka a 

+ ^Y,nln f _ a + Vj2(cU« + rtck a ) (21) 

a ka 

consists of a single local state, which we will denote with 
/, with energy e/ and Coulomb repulsion U, coupled to 
a bath of conduction electrons with creation operators 
c\ a and energies Ck a - The local level is subject to a Zee- 
mann splitting in an external magnetic field B. Note 
that we denote with H = gfi B B/2 the Zeeman energy 
and the total splitting | e-p — ej = 2H. We employ Wil- 
son' s NRGi 2 ^ to generate the eigen-energies, the matrix 
elements and basis set needed for the Green's function. 

So far, all analytical calculations were performed using 
the discrete NRG spectrum. To obtain a continuous spec- 
tral function from the set of discrete ^-functions occurring 
in Ga,b(z), we have to introduce a coarse-graining. Due 
to the exponentially decreasing energy scales, a broaden- 
ing on a logarithmic mesh by a Gaussian 

is typically usedii^S Since this broadening function is 
properly normalized to one, the spectral weight is con- 
served and no principle inaccuracies are introduced by 
this procedure. 



B. Spectral functions for T = 

In the following, we specialize to A = / CT and B = /J. 
For comparison we calculated the raw NRG spectral func- 
tion in three different ways: (i) by the conventional way 
as described in Refs. IT5IIT7I labelled as (CON) in the fol- 
lowing, (ii) by the Hofstetter approach^ (HA) and by 
our method defined by Eqs. (|ll|) - (|17|l labelled complete 
Fock space approach (CFS). As long as not stated other- 
wise, all energies are measured in units of the half band 



width D, and for simplicity only a symmetric conduc- 
tion band is considered with a constant density of states 
po = l/(2£>)0(£>-H)Jl 
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FIG. 2: (color online) Comparison of the spectral function for 
the three different methods for the symmetric case tf/D = 
— U/D2 = 0.5. The inset on the left side emphasizes the broad 
charge excitation at £/, the inset on the right side zooms in 
at the Abrikosov-Suhl resonance (ASR). NRG parameters: 
F/D = irV 2 p /D = 0.1, A = 2.5, N s = 500, b = 0.8, T = 

In Fig |2 we show a comparison of spectral functions 
calculated by the three different methods using identical 
NRG input data and broadening parameters. All three 
methods agree on the shape and the height of the many- 
body resonance (ASR) at the Fermi energy and all meth- 
ods do not reach the unitary limit of TrTp(0) = 1 as pre- 
dicted by the density of states rulei2£ The spectral func- 
tions, however, differ in their high energy features. The 
conventional spectral function underestimates the charge 
excitations and only reaches a total spectral weight of 
C = 0.946. The Hofstetter approach^ yields a somewhat 
improved value of C = 0.956 because the high energy 
part of its spectrum, which carries most of the spectral 
weight, lies between our and the conventional curve. The 
spectral function based upon our new approach, however, 
reaches the spectral weight C = 1.00034. It deviates from 
the exact value of C = 1 only by the error introduced by 
the numerical uj integration. To verify this, we add all 
spectral weights from Eqs. I|ll|) - (|16[l directly and obtain 
\C - 1| = 10" 12 . . . 10~ 15 , i.e. the norm is C = 1 within 
machine accuracy42. For this symmetric case, its obvious 
that our approach also yields the exact occupation num- 
ber of nl = (hi) =0.5 for each spin direction while the 
other two approaches have a 5.7% (CON) or 4.6% (HA) 
error. Note, however, that the latter two still give the 
correct value C/2 with respect to their respective norm. 

At this point a comment on the violation of the density 
of states rule 7rFp(0) = 1 for T — > seems in order. Since 
the NRG itself does only provide the weight of <5-peaks, 
the height of the coarse-grained spectrum will depend 
on the density of these peaks, i.e. the number of states 
available in the calculation, and the choosen broadening 
function. This does not imply the violation of Friedel 
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sum rilled 9 ., too, which relates the scattering phases to 
the number of displaced electrons in each spin channel. 
Since the NRG fixed point spectrum contains the correct 
scattering phases, it is accurately fulfilled. However, the 
relation between the phase shift and the local spectral 
function is in general not very accurately reproduced due 
to the aforementioned reasons. Nevertheless, the density 
of states sum rule can be satisfied with an error of less 
than 3% using small values of A, more NRG states and 
a smaller broadening^ as used in Fig [21 

To our knowledge, only the indirect calculation of the 
physical Green's function by expressing its self-energy 
by the ratio of two correlation functions^ yields val- 
ues which are accurate enough to reproduce the correct 
height of p(0, T = 0) nearly independent of the broad- 
ening parameters. In particular, this is of importance to 
prove the pinning of the spectral function for the two- 
channel Anderson model at half the unitary limit. 23 
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FIG. 3: (color online) Comparison of the spectral function 
for the three different methods for the same model and NRG 
parameters as in Fig. [5] but with only Ns ~ 50 states kept in 
each NRG iteration. The inset zooms in onto the ASR. 

We have put our method to an extreme test by reduc- 
ing the number of states Ns kept after each NRG step 
to Ns = 50 for the same A = 2.5 as before. This is def- 
initely stretching the limit of the conventional NRG be- 
yond its accuracy. The results for the spectral functions 
are plotted in Fig.|3] Only data for frequency |u/| > 10~ 7 
are included, which, however, is much smaller than the 
Kondo temperature Tk as estimated from the width of 
the ASR. The inaccuracy of the matrix elements for ex- 
tremely small frequencies \u>\ < 10 ~ 7 yields erratic re- 
sults in this interval which, therefore, are omitted here. 
While the conventional and the HA methods significantly 
loose spectral weight, the spectral weight of CFS spec- 
tral function remains at the exact value C — 1. Since 
the discretization of the energy mesh was the same in all 
figures, the deviation from the sum rule of 3 x 10~ 4 again 
stems only from the numerical integration. For compari- 
son we also added the CFS curve from Fig. [2 calculated 
with N s = 500 NRG states. Both CFS curves agree ex- 
tremely well. In particular, the weight under the charge 



peak is only slightly redistributed and the shape of the 
ASR remains unaltered. Only at very low frequency, the 
inaccuracy of the underlying NRG input data takes its 
toll: the ASR peak height is further reduced. The other 
two methods, however, exhibit strong dependence on Ns 
in all energy regimes. We must report an error in the 
spectral weight and the occupancy w a of w 35% for both 
the conventional and the HA method. In all cases the 
same broadening parameter b = 0.8 was used. 




CO 



FIG. 4: (color online) Comparison of the spectral function 
obtained from CFS and HA approaches in a finite magnetic 
field (a) H = 0.005 and (b) H = 0.1. The full curves represent 
spin up, the dashed ones spin down. The inset in (a) is a zoom 
into the region around the Fermi level showing the splitting 
and suppression of the Kondo resonance. Model and NRG 
parameters as in Fig. [3] The norm and occupancies obtained 
are compared in Tab. [I] 

While the conventional method yields quite reliable 
spectral information for small A and large number of 
states in the absence of a magnetic field, it fails for finite 
magnetic field, where the fixed point for T — > yields a 
complete redistribution of spectral weight on all energy 
scales compared to the spectral function^ at H = 0. 
Therefore, it is omitted in Fig.^ where we compare spec- 
tral functions for two magnetic fields, H = 0.005 s» Tk 
(Fig. Hi) and H = 0.1 > T K (Fig. @J>); the results for 
the majority spin are plotted as solid lines, the ones for 
the minority spin as dashed lines. For H w Tk the ASR 
is split and already reduced to half its original height, 
while for H ^> Tk it has vanished completely, as ex- 
pected. Again, we display the data for Ns — 50 to high- 
light the differences in the methods. While our method 
and the HA approach exhibit the same overall features, 
i.e. splitting and reduction respectively complete lack of 
the ASR and a major redistribution of spectral weight, 
the charge excitation peak contains much more spectral 
weight in the CFS as the HA curve. As before, the CFS 
spectral function fulfills the spectral sum rule with the 
accuracy of 10 -4 while we note a 20% error in the HA 
curve for these extreme set of parameters. The CFS re- 
sults for Ns — 500 given by the blue curves in Fig. 0] by 
and large again lie on top of the CFS curves for Ns = 50, 
once more emphasizing that within the CFS neither the 
accuracy of the occupancy nor the one of the spectral 
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Method 


Norm 




)gf 


(net) 


NRG 








a=| 


0=1 


a=| 


0=1 


(a) 


CFS(500) 


1 


0.807 


0.193 


0.807 


0.193 


H = 


0.005 CFS(50) 


1 


0.806 


0.194 


0.806 


0.194 




HA 


0.747 


0.600 


0.147 


0.806 


0.194 


(b) 


CFS(500) 


1 


0.953 


0.047 


0.953 


0.047 


H = 


0.1 CFS(50) 


1 


0.953 


0.047 


0.953 


0.047 




HA 


0.793 


0.755 


0.038 


0.953 


0.047 



TABLE I: Norm and occupancies for the calculations with 
finite magnetic field in Fig. |1] The last columns contain the 
occupancies obtained from the thermodynamic expectation 
values. 



functions does critically depend on the number of states 
kept. 

This feature becomes even more apparent from the ac- 
tual numbers for the norm and the spin-dependent oc- 
cupation numbers listed in Tab. HJ The last column 
shows for comparison the occupancies for the two spin 
directions as calculated directly from the NRG level 
flowiii Note that the difference in the occupancies be- 
tween Ns = 500 and Ns = 50 in the CFS appear in the 
thermodynamic occupation numbers, too. On the other 
hand, there exist significant differences in the magneti- 
zation m = (ri| — ni)/2 obtained from the HA method, 
which yields a 30% error compared to the reference NRG 
magnetization and our CFS method. 

We like to emphasize that the rather large error of the 
HA method results from the unusual low number of NRG 
state Ns kept after each NRG iteration. Increasing the 
number of states and reducing A one can increase the 
accuracy of the HA method considerably, but typically 
is always left with about 4% error according to Tab. I in 
Ref. UJJ. Our motivation for choosing Ns = 50 was to put 
our initial claim to an extreme test: The spectral sum 
rule and the occupation numbers are reproduced with 
high accuracy independent of the number Ns of states 
kept after each iteration. Our equations \ll\) - \lb)) do not 
contain any truncation errors since we use a complete 
basis set. The errors in the CFS spectral functions are of 
purely energetic nature and will indeed be found in the 
thermodynamic quantities, too. As a consequence, the 
spectral functions appear to be largely insensitive to the 
number of NRG state kept. 

We made the same observation as discussed above 
in the asymmetric case shown in Fig. Efa). The over- 
all shape of the spectral functions of all three methods 
agree pretty well. Again, we note differences in the de- 
viations from the sum rules, see Tab. [H] While our 
new CFS method obeys the spectral weight sum rules 
and the occupation sum rule very accurately, the com- 
parison of the actual numbers for norm and occupancy 
in Tab. |Hja) shows that the conventional and the HA 
method yields C — 0.974 and the Green's function occu- 
pancy of n^ F = 0.645 deviates from the NRG occupancy 
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FIG. 5: (color online) Comparison of the spectral function 
for the three different methods for the asymmetric case (a) 
U/D = 1 and e f /D = -0.9, and (b) U/D = 1000 and 
ej/D = -0.4. NRG parameters: F/D = TvV 2 p /D = 
0.1, A = 2.5, Ns = 500. The results for norm and occupancies 
obtained are collected in Tab. [H] 





Method 


Norm 


(n)oF 


{n) NRG 


(a) 


CFS 


1 


0.666 


0.666 


U = 1 


HA 


0.974 


0.645 


0.666 




CON 


0.975 


0.645 


0.666 


(b) 


CFS 


0.543 


0.915 


0.915 


U = oo 


HA 


0.523 


0.878 


0.915 




CON 


0.534 


0.893 


0.915 



,NRG 



0.666 by 3%. 



TABLE II: Norm and occupancies obtained with the different 
methods to calculate spectra for the asymmetric SIAM (see 
Fig. for the parameters used). The last column shows the 
occupancy as obtained from the thermodynamic expectation 
value. 



In Fig. Etb) the spectral functions for U/D = 1000 are 
plotted. This corresponds to the limit U — > oo of the 
SIAM, i.e. the charge excitation between the singly and 
the doubly occupied state is effectively shifted to infin- 
ity. Therefore, the total spectral weight is reduced to 
C = 1 — (hf)/2, if we neglect the <5-like peak at to w U. 
Nevertheless, adding all discrete NRG spectral weights 
yields |C — 1| ~ 10~ 12 — 10~ 15 even in this case. Again, 
the values in Tab.lTlTb) show that the CFS spectral func- 
tion agrees excellently with the exact norm but show a 
slight error of 0.03% in the occupancy, which we again 
attribute to the numerical integration, while the conven- 
tional and the HA methods have errors at least two orders 
of magnitude larger of about 4%. 



C. Finite temperature spectral functions 

Let us now turn to the discussion of the temperature 
dependence of the spectral functions calculated with CFS 
method. Here, an additional problem arises, because the 
temperature scale in NRG is defined by the length N of 
the chain. More precisely, a given chain length N cor- 
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FIG. 6: (color online) Spectral function for temperatures T = 
0, 0.00453 w Ik, 0.01133 and 0.112 > T K calculated with 
CFS. The full curves result from calculations with Ns = 500, 
the dashed ones from calculations with Ns = 50. Model and 
NRG parameters are otherwise the same as in Fig. 



responds to a temperature Tn ~ f^-( N - 1 )/' 2 _ Thus, a 
continuous variation of T is not possible, only a discrete 
set of logarithmically varying temperatures is accessible. 
More problematic, however, is that due to the termina- 
tion of the NRG iterations at chain length N no infor- 
mation are available for excitation energies smaller than 
TV, too>i£ This means, that the part of the spectrum 
with \uj\ < Tjv is inaccurate in the sense that excitations 
on these scales cannot be taken into account properly. 
Furthermore, the limited information available for this 
energy range stems from all NRG iterations, and hence 
a Gaussian broadening l|22fl suitable for excitations col- 
lected from a single iteration cannot be used here. In- 
stead, a Lorcntzian broadening 15 



1 



27T ( w - Wn )2 + b 2 



(23) 



is used for |cj| < aT/v, where a sets the energy scale 
down to which we trust the NRG results^ Moreover, as 
discussed in detail in the appendix, the CFS still ensures 
a complete basis set and hence an exact norm for the 
spectral functions, but quantities obtained from spectral 
averages like occupancies must be expected to be less 
accurate than for T = 0. 

The development of the spectral function for T/D = 
ttV 2 p/D = 0.1, e f = -(7/2 = -Q.5D as function of 
temperature can be found in Fig. Even though we re- 
strict the results presented in this section to the particle- 
hole symmetric limit, all observations also hold for the 
asymmetric case with and without magnetic field. As 
expected, with increasing temperature the ASR is re- 
duced, for T sa Tk to roughly half its original height. 
For T ^> Tk it has completely vanished, leaving only 
the Hubbard peaks in the spectrum. For all spectra we 
find the norm to be exactly one as before and, due to 
particle-hole symmetry, the occupancy has to be one- 
half. The full curves in Fig. are results of NRG cal- 



culations with Ns — 500 states, the dashed ones with 
Ns = 50 states. As before, we can notice mild devia- 
tions of the spectra in the region of the Hubbard peaks 
and at very low energies, but the overall agreement of the 
spectra is quite good, thus again underlining the previous 
claim that spectra calculated with the CFS are rather in- 
sensitive to the number Ns of states kept in each NRG 
iteration. Note that the deviations for u> — > can be eas- 
ily accounted for by the fact that with varying Ns the 
distribution of excitation energies in the NRG will differ, 
thus yielding a slightly different distribution of spectral 
weight due to the broadening lj2l5|) . 

The calculations for finite magnetic field H — 0.005 
are collected in Fig. [7| Here we only present results for 



— T=0.112 

— T=0.01133 

— T=0.00453 

— T=0 




FIG. 7: (color online) Spectral function for temperatures T = 
0, 0.00453 » T K , 0.01133 and 0.112 > T K and finite magnetic 
field H = 0.005 calculated with CFS. The full curves show the 
spectral function for spin up, the dashed curves spin down. 
Model and NRG parameters are otherwise the same as in 
Fig. 1^1 The values for norm and occupancies are listed in 
Tab. EH 

Ns — 50, because as noted before the spectra for Ns = 
500 do not differ significantly. As in Fig. [3] we plotted the 
spectra for the majority spin as full curves, the ones for 
the minority spin as dashed curves. Again, the norm is 
exactly one in all cases. The occupancies obtained from 
integrating the spectra multiplied with the Fermi func- 
tion at the appropriate temperature are listed in Tab. lllll 
for Ns = 500 and Ns = 50. To avoid additional inaccu- 
racies introduced by the broadening (|23[) we have listed 
in Tab. IHII the sum of the weights of the raw NRG spec- 
tra multiplied with the Fermi function. As for T = 0, 
the deviations to the thermodynamic values stay always 
less than 0.01%. Numerical integration of the broad- 
ened spectra, on the other hand, will lead to stronger 
deviations, in particular for higher temperatures. For 
example, we find (n{)cF — 0.516 for T = 0.112 and 
Ns = 500, i.e. an error in the percent range. The value 
of J p{ui)f{uj) with p(ui) smoothened by the described 
combination of Gaussian and Lorentzian broadening de- 
viates from the sum over the spectral weights of the poles 
of the discrete spectrum times the Fermi function at the 



10 



Ns 


T 


Norm 


(n a 


)gf 


(no-) 


NRG 








a=T 


0=1 


a=T 




500 





1 


0.806 


0.194 


0.806 


0.194 


50 





1 


0.807 


0.193 


0.807 


0.193 


500 


0.00453 


1 


0.723 


0.277 


0.723 


0.277 


50 


0.00453 


1 


0.729 


0.271 


0.729 


0.271 


500 


0.01133 


1 


0.629 


0.371 


0.629 


0.371 


50 


0.01133 


1 


0.635 


0.365 


0.635 


0.365 


500 


0.112 


1 


0.518 


0.482 


0.518 


0.482 


50 


0.112 


1 


0.519 


0.481 


0.519 


0.481 



TABLE III: Norm and occupancies for the calculations with 
finite magnetic field and T > in Fig. [7| In addition the 
occupancies for Ns = 500 are included in the table. The last 
columns contain the occupancies obtained from the thermo- 
dynamic expectation values. 



pole energy, since the Lorentzian part extends in an un- 
controlled manner to positive and negative frequencies. 
This systematic error, however, is still small even in these 
extreme cases. 



IV. CONCLUSION 

We presented a novel approach for the calculation of 
local Green's functions for quantum impurity systems us- 
ing Wilson's numerical renormalization group. Since it 
is based on the usage of a complete basis set for the Wil- 
son NRG chain, recently introduced in the context of the 
time-dependent NRGj22*2I quantities written as spectral 
sum rules are always exactly fulfilled, regardless of the 
size Ns of the Hilbert space kept in each NRG iteration. 
Moreover, spectral averages like occupation numbers etc. 
are at least reproduced with unprecedented accuracy. We 
have demonstrated the quality of the approach by pre- 
senting results for the fermionic single-particle spectral 
function pf{uf) of the single impurity Anderson model. 
We have shown explicitly the validity of our claim, that 
the resulting spectral functions are very insensitive to the 
number of NRG states kept in the iteration. This is a con- 
sequence of using a complete basis set in the derivation 
of the Lehmann representation. Equations l|lll) - (|16fl P rc ~ 
cisely account for which excitation contributes at which 
energy shell m, thus manifestly solving the double count- 
ing problem of excitations responsible for the violation 
of the sum rules in all previous NRG spectral function 
approaches: The conventionally employed "patching" al- 
gorithms thus have become obsolete. Also, no even-odd 
oscillations can be observed. We believe that our novel 
approach provides the most accurate representation of 
local spectral functions for a given NRG input. It can be 
used in quantum impurity systems, at finite temperature 
as well as T — > 0. It is particular suitable for symmetry 
broken phases, for instance in an external magnetic field. 

This significant improvement can have a major impact 
on the usage of the NRG as impurity solver for the dy- 



namical mean field theory (DMFT) AS In particular, for 
the description of symmetry broken phases such as the 
ferromagnetic Hubbard models as well as orbital order in 
the two-band Hubbard model, the discrepancy between 
the local NRG order parameter and the one obtained 
from the spectral function often yields instabilities in it- 
erating the DMFT equations. Moreover, very accurate 
spectral functions are needed for the application of the 
NRG to cluster DMFT approaches^ or periodic Ander- 
son models which include crystal field effects^. Up to 
now such an accuracy could only be achieved with im- 
mense CPU time4£ 
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APPENDIX A: PROOF OF THE SPECTRAL 
SUM RULE 

By inserting the analytical expression (|17|l for the dis- 
crete spectral function into the integral (|19fl 

C = /S G "M= £ /s^w (A1) 

yields the following three contributions 

C = J2 A H'( N ) B i'A N ) (e- pE ' -se-f"®) (A2) 
l,l< 

N-l 

+ E EE 5 ^'( m )^( m ) A ^w 

m—m m i n I k.k' 
N-l 

- s E EE^w^(*h • 

m=m m in I k,k' 
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By formally defining the "reduced density matrix" for 
the last Wilson shell as 



pltf(N) = 5 k , k , 



(A3) 



the first term in l)A2jl can be included into the second 
and third. Then the second reads 



/v 



m— m„ 
N 



I' k,k' 



= E EEE< r ' e ";™w,e';m> 

m=m min I'e" k,e k' ,e' 

x (fc', e'; m\p\k, e; m)(k, e; to|A|Z', e ; m) (A4) 

In this step, we used the definition of the reduced density 
matrix and made use of the fact that the local operators 
A and B are independent of the environment variables 
e, e', e". Due to the form of the density operator, (1+ + 

lm)p(lm + lm) = im^m holds for 171 < N - For A = TO, 

ljy spans the complete Fock space and the projection 
l~pl~ contains all contributions at m = N. Therefore, 
the contribution C2 yields 



JV 



C 2 = ]T ^{l',e";m\B P A\l',e";r 

= Tr [pAB] 



m=rn min I'e" 



(A5) 



We can perform the same calculation for the third term 
in (|A2(I and the contribution stemming from the (— s) 
part of the first term at m = JV to derive 

C* 3 = -sTr [pBA] . (A6) 

To this end, the total contribution is given by 

C = Tr [pAB] - sTr [pBA] = Tr [p[A, B] s ] (A7) 

We thus have proven that our Green's function fulfills the 
spectral sum rule exactly. Thus, any deviations for this 
sum rule in the broaden spectral function stems solemnly 
from the accuracy of the numerical oj integration. 

APPENDIX B: OCCUPATION SUM RULE 

We now want to discuss the accuracy of our approach 
for the generalized occupancy ub.a defined as 



ns.A 



dz 
2iri 



f s {z)G A<B {z) 



(Bl) 



where f s {z) = [exp(/3z) — s] _1 . As always through- 
out the paper, this definition comprises bosonic and 
fcrmionic Green's functions. Specializing to the local 
spin-dependent fermionic spectral function, i.e. A = f„ 
and B = f£, yields the occupational sum rule i|2(JI) as 
presented in section III Dl 



Inserting the Green's function l|17|l into the expression 
(|B1() and performing the contour integration yields the 
following three contributions: 



-PEl 



riB, 



a = ^2^-Bi,i>{N)Ai,,i(N) 



N-l 



(B2) 



2^ 2^2^, e 0(E r - E ^) _ s 

m—nimin I k,k' 



N-l 



\ - \ y A hk ,(m)pf, d k (m,)B k: i(m) 

2^ 2^2^ l _ se p( E ™-E™) 

m=m mirl I fc,fc' 



Adding and subtracting 



N-l 



E EE 

m—rrirnin I k.k' 



Ak'(m)p k e , d k (m)B kil (m) 



(B3) 



leads to 

n B ,A = ^2—^-B, tV (N)Ai f ,i(N) 
l,l' 

N-l 

+ E EE^^'H^h^'H 

m—mjnin I k,k' 

+Stib,a 
= Tr [pBA] + 5n B ,A , 

using the same identities as in appendix Finally, the 
error Sub, A of the occupation sum rule is given by 



Sub, 



N-l 

E 

m—rrimin I, k.k 



E P^%(m) (B4) 

,k,k' 

B Lk '(m)A kt i(m) - A iik > (m)B k ^(m) 
X e ^ E r- E T) - s 

To estimate the size of this error, let us note that the exci- 
tation energies EJ 71 — E™ appearing in Sub, a are positive, 
because I labels a discarded and k a kept state. Excita- 
tions between two discarded states never contribute as 
well as excitations between to kept states. The small- 
est energy difference possible is given by the energy of 
the first discarded minus the energy of the last kept 
state. Even though this energy difference might be 
small, the smallness of the matrix element of the re- 
duced density matrix suppresses this contribution. We 
are thus left with energy differences which are of the or- 
der D m for to < A. The density matrix p, on the other 
hand, is restricted to the last iteration, i.e., to be con- 
sistent with the fundamental assumption of the NRG, 
exp[/3(£; ; m - Eg)] > 1 for to < JV. Thus, the denom- 
inator in 8riB,A becomes exponentially large and hence 
Sns, a exponentially suppressed with decreasing temper- 
ature. Note that this also means, that for T — > the 
occupation sum rule will be obeyed exactly, too. For fi- 
nite T we must expect, and indeed find (see section lTlI CL 



12 



deviations from this sum rule with increasing tempera- 
ture. Our numerical results show, however, that even for 
high temperatures these deviations stay below the per 
mille level. 



APPENDIX C: IMAGINARY TIME GREEN'S 
FUNCTIONS 

The textbook definition of the imaginary time Green's 
function is given by 

GaM-t) = -Tr [pT (A(t)B)} 

= { -Tr [pe TH Ae- TH B] ; r>0 
[ -sTr [pBe rH Ae~ TH ] ; r < 

where again s = 1 for Bosonic operators A, B and s — 
— 1 for Fermionic operators. It is straight forward to 
show that the analytic continuation Ga,b{z) of the exact 
Matsubara Green's function Ga,b{^) 

Ga,b(^) = / dTe^ T G A , B (T) , , (C2) 
Jo 

where to = n(2n + l)//3 for Fermionic and u> = 2i:n/(3 
for Bosonic Green's functions, is identical to the Laplace 



transformed exact retarded Green's function with z = 
iuj. We can perform the same steps for Ga,b(t) as in 
section HTGl In this case, however, the terms — s[z + 
E m _ E™Y X in G u {z) and [z + Eg - E™Y X in G m {z) 
must be replaced by 



exp(/?(£™ - Ef 1 ) - s 
iu + E™ - E™ 



and 



iuj + E™- E™ 

At first sight, the two Green's functions are not equal in 
our approach. However, as discussed in appendix FBI the 
Boltzmann factors exp[/3(£'™ — E™)] can be neglected for 
m < N. This is justified by the NRG assumption that 
the density operator is well approximated by its contri- 
bution of the last Wilson shell, i.e. Eq. 0). Therefore, 
both Green's functions are identical within the NRG ap- 
proximation. 
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We use the standard notation, by which the AT-site chain 
contains the impurity degrees of freedom, as well as the 
first N + 1 Wilson shells (labelled by n = 0, ■ ■ • ,N) 
In contrast to the standard NRG notations we consider TLn 
as dimensionfull Hamiltonian in order to kept the notation 
as simple as possible. 

Local operators we call operators which act only on chain 
links up to m < m m i n . See a detailed discussion in Refs. 
I2(3l2l 

We can either directly sum the weights of the discrete raw 
NRG spectrum or perform a numerical integration after 
having broadened it. The former usually yields an order of 
magnitude better accuracy than the latter. 
This in fact holds for all spectra calculated, i.e. with or 
without magnetic field, at T = and finite temperature! 
For our calculations, we choose b = 0.8 and a — 2. 



